Comments on Bona-Masso type slicing conditions in long-term black hole evolutions 
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We review in generality why time-independent endstates can be reached in black hole and collapse 
simulations, with and without excision. We characterise the Killing states of the Bona-Masso slicing 
condition with time derivative along the normals to the slice ("BMn") as solutions of a mixed 
elliptic/hyperbolic differential equation on the slice. We show numerically that these steady states 
can be reached as end states from typical initial data with excision but can be reached with the 
puncture method only if the puncture is not numerically well resolved. During the evolution, BMn 
slicings often form gauge shocks. It may be that these are not seen in current 3D simulations only 
through lack of resolution, although we expect that they can be avoided with some care. Finally we 
point out that excision with BMn as currently implemented is ill-posed and therefore not expected 
to converge; this can be cured. In technical appendixes, we derive the equations of pure gauge 
systems on a fixed spacetime, and bring the BSSN/NOR equations into 3-dimensional tensor form 
suitable for multiple coordinate patches or spherical polar coordinates. 
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I. INTRODUCTION 

In 2005, parallel breakthroughs in the long-term sta- 
ble simulation of binary black holes were made using 
two rather different approaches: Pretorius [l| used modi- 
fied harmonic slicing, black holes created in collapse, and 
singularity excision, while the Brownsville '^l and God- 
dard [3] groups used eternal black holes, and singularity- 
avoiding (l-|-log) slicing. Nevertheless, a key ingredient 
in both these successes is their gauge choice. 

Generalising and extending recent work by Hannam et 
al. [4|, we further investigate the application to black 
hole and collapse spacetimes of the Bona-Masso slicing 
condition with time derivative along the slice normals 
("BMn"). This family includes both the "1-Mog" shcing 
used in [3j3j| and the harmonic slicing, a variant of which 
is use in 

A desirable property for a gauge choices is that the 
metric becomes time-independent to the extent that the 
spacetime becomes stationary Q. In Sec. |n]we explain 
carefully why this is possible both when the black holes 
are excised and when a singularity-avoiding slicing is 
used. We characterise Killing coordinates geometrically 
in Sec. IIIII To fix notation, we review various lapse con- 
ditions of Bona-Masso type in Sec. lIVI In Sec.|V]we clas- 
sify Killing slicings compatible with BMn slicing, and in 
particular the spherical Killing slicings of Schwarzschild 
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spacetime. In Sec. IVII we investigate numerically if any 
such Killing states are in fact attractors in evolutions 
of the Schwarzschild spacetime. We consider both slices 
with wormhole topology and slices which end at an ex- 
cision boundary inside the black hole. In Sec. IVIIl we 
present spherically symmetric simulations of scalar field 
collapse as a toy model for black holes formed in collapse. 
From our mathematical and numerical observations in 
these sections, we suggest improvements to current meth- 
ods for binary black hole evolutions in Sec. IVIIII 




II. NUMERICAL EVOLUTION OF BLACK 
HOLE SPACETIMES 

A. Eternal and collapse black holes 

Black holes in the real world have formed in collapse, 
but eternal black holes are often used in numerical rela- 
tivity because they differ from collapse black holes only 
in the interior, and this cannot affect physics outside. 
Here we concentrate on non-rotating, uncharged black 
holes, which are described by the Kruskal extension of 
the Schwarzschild spacetime. A bifurcate Killing horizon 
divides this spacetime into past (P), future (F), "left" (L) 
and "right" (R) regions. The future and past timelike 
and i~) and null and J^~) infinities and the space- 
like infinity i° all exist in left (L) and right (R) copies. 
Slices extending from j° to have wormhole geometry, 
see Fig. [TJ Binary (or multiple) black hole initial data 
can be represented by a wormhole leading to a separate 
copy of for each black hole. In the "puncture" method 
, each is then represented in coordinates by a point 
where the conformal factor diverges. 

By contrast, black holes formed from regular data 
through collapse have trivial spatial topology, similar to 
the Schwarzschild spacetime but with part of R and F, 
and all of P and L, covered up by the collapsing star Q 
- see Fig. |4l 

B. Singularity-avoiding slicings 

Both in collapse and in eternal black holes one can 
use slicings which avoid the singularity. Any timelike 
worldline inside a black hole has finite length, while any 
timelike worldline with limited total acceleration outside 
the black hole has infinite length. The lapse measures the 
rate of proper time per coordinate time for an observer 
normal to the time slices, and so one might think that 
the lapse must go to zero everywhere inside the black 
hole in order to avoid the singularity, and that because 
the slices keep advancing outside the black hole, their in- 
trinsic geometry must deform without limit as time goes 
on, until large gradients can no longer be resolved. Such 
"slice stretching" was indeed encountered in early black 
hole simulations, and motivated the development of black 
hole excision [8i]. 



FIG. 1: Spacetime diagram of the Schwarzscliild spacetime, 
with the angular coordinates suppressed. The horizontal line 
from i° to i% is the time-symmetric wormhole slice typically 
used as initial data in puncture evolutions of a Schwarzschild 
black holes. The curved lines schematically represent the slic- 
ing generated from these initial data by BMn lapse with a = 1 
initially. They approach the slice R = Rq, which links to 
ij. The vertical dashed line represents the symmetry bound- 
ary which can replace the left-right reflection symmetry of 
this slicing. As the slices approach R — Ro, the approxi- 
mately cylindrical wormhole grows longer linearly with time. 




FIG. 2: The same spacetime diagram, schematically showing 
the unique regular spherical Killing slicing that is compatible 
with BMn slicing (for a given /i_t(Q). All slices are isometric 
to one another, and connect with t^. The again asymptote 
to the slice R = Rq. 

Only later it was realised clearly that singularity- 
avoiding slicings need not lead to slice stretching 
If the lapse is chosen such that the slice is Lie-dragged 
along the Killing vector field everywhere, its intrinsic ge- 
ometry becomes time-independent. This is true also in- 
side the black hole where the Killing vector field that 
generates time translations at infinity becomes space- 
like (and so the spacetime is not technically station- 
ary), as long as this Killing vector field is nowhere par- 
allel to the slicing. Once the geometry of the slice 
has become time-independent, a suitable shift condition 
then makes the spatial metric coefficients explicitly time- 
independent. With this lapse and shift d/dt becomes the 
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FIG. 3: The same spacetime diagram, schematically showing 
a Killing slicing that ends at the future singularity, such as 
Kerr-Schild slices. The lines with arrows are trajectories of 
the Killing vector (lines of constant R) and the beads on them 
represent surfaces of constant coordinate r if the Killing shift 
is used. In particular, the dashed line could serve as a Killing 
excision boundary. 




FIG. 4: Schematic spacetime diagram of the collapse of a 
spherical star. Outside the collapsing star (shaded) the space- 
time is Schwarzschild, comprising parts of regions R and F. A 
Killing slicing with excision as in Fig. [3] is shown. A Killing 
endstate cannot be reached without excision. 



Killing vector (spacelike inside a black hole) . Coordinate 
conditions which generate Killing coordinates asymptoti- 
cally starting from generic initial coordinates were called 
"symmetry-seeking" in [5|. 

Even more recently it was realised that the lapse need 
not collapse either Q. Note that 



is a sum of two terms. Define some scalar a to measure 
distance from the singularity. (In Schwarzschild space- 
time, an obvious choice is the area radius R.) For any 
given a and cr, /3* then can be chosen to set a = 0, ex- 
cept where cr ,; = 0. (We use a dot to denote d/dt). 
In other words, the lapse in a Killing coordinate system 
vanishes only where the time slices are tangential to the 
Killing vector field. Every regular time slice in a collapse 
spacetime, and every wormhole slice through an eternal 
black hole has such an obstruction point, namely a local 
minimum of a (Fig. [1]). However, a slice that becomes 
asymptotically cylindrical (with R ^ Rq) and ends at 
avoids this obstruction (Figl2]). 

C. Excision 

An alternative to singularity-avoiding slicings is sin- 
gularity excision. This means truncating the time slices 
along a future spacelike surface which is also (at least 
asymptotically) Killing. In Schwarzschild spacetime, this 
would be a surface of constant R < 2M. One still wants 
the slice to be Lie-dragged along the Killing field, but 
one gains more freedom because Killing slices are now 
acceptable which would intersect the singularity, such as 
Kerr-Schild slices of Schwarzschild. A Killing slicing with 
Killing excision boundary is illustrated in Fig. [3l 

As long as the excision surface is spacelike, all char- 
acteristics corresponding to gravitational waves, which 
propagate on light cones, will be leaving the domain of 
computation. Depending on the formulation of the Ein- 
stein equations and the gauge choice, other characteris- 
tics corresponding to constraint modes and gauge modes 
may be spacelike, and either this will restrict the exci- 
sion surface further or explicit boundary conditions need 
to be imposed on the gauge if the evolution equations 
are to be well-posed. If the system is not hyperbolic, 
for example because the gauge conditions are parabolic 
or elliptic, boundary conditions will be required on any 
excision bovmdary. 

In [10] and JTj , evolutions were carried out from punc- 
ture data using BMt 1-1- log slicing with and without Kq, 
and directly comparing evolutions using either excision or 
fixed punctures. No explicit boundary condition was im- 
posed at the excision boundary. Excised and non-excised 
evolutions are claimed to converge to each other to sec- 
ond order everywhere outside the excised region. This is 
surprising given that the excision problem was ill-posed. 

III. KILLING COORDINATES 

A. General case 

By definition, coordinates in which the 4-metric is 
time-independent are those in which {d/dtY — C*^", 
where is a Killing vector that is timelike at infinity 
and C 7^ is a constant. Contracting with ria, we find 
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that the Killing lapse is given by 

a = C<j>, (2) 

where (j) = —na£,°', and contracting with the projector 
_L(j^ = Qa'^ + n^'n!' we find that the Killing shift is 

/3' = ?aer- (3) 



B. Schwarzschild spacetime in spherical symmetry 

We now restrict to spherically symmetric Killing coor- 
dinate systems on the Kruskal extension of Schwarzschild 
spacetime. In the following, are preferred coordinates 
on a given spacetime such as Schwarzschild, while it, x^) 
are the coordinates used for the numerical evolution, in 
our case with the spherical line element 

ds^ ^-a^dt^ +-^{dr + (3dtf +R^dn'^, (4) 

We use the shorthands dVl^ = df?^ + sin^ 9dip^, R"^ = ^gg, 
7 = ^rr and /3 = /J''. We use / and /' for the partial 
derivatives with respect to t and r. 

We use preferred coordinates (T, R) on Schwarzschild 
with the property that R is the area radius and the 
Killing vector is d/dT, normalised to unity at infin- 
ity, for example Schwarzschild or Kerr-Schild coordi- 
nates. In all such coordinates grr = 1 — 2M/R and 
grrgRR. — gxR = — l-The generic Killing coordinate sys- 
tem {t, r) with C = 1 is then given by the ansatz 

T = t + F{r), R = R{r). (5) 

If we are interested only in the slicing, we can fix the 
spatial coordinate r for convenience. A better choice than 
using R itself as a coordinate is to make r proper distance 
along the slice, so that 7 = 1. (We shall also use the 
symbol I for proper radial distance.) The Killing lapse 
and shift are 

a = R', (6) 

(3 = + (7) 

The trace of the extrinsic curvature of the Killing slices 
is 

13' 

+ (8) 

where /3 is given by Q. 

IV. EVOLVED SLICING CONDITIONS 

We focus on the family of slicing conditions suggested 
by Bona and Masso [12] (from now BM) 

an°-V aa = a - (3'a,i ^ -^iLa^K, (9) 



where K is the trace of the extrinsic curvature of the 
slice and its unit normal vector. Typically, /ii, > is 
understood to be a given function ^L{ct) of the lapse. As 
n"" is a true vector and a and K are scalars under a change 
of coordinates a;* on the slice, this slicing condition is 
independent of the coordinates on the slice and therefore 
independent of the shift. 

Confusingly, the very different slicing condition 

a^-tiLa^{K~Ko), (10) 

where Ko{x^) is the initial value of K [9|], is also referred 
to as Bona-Masso slicing. For clarity, we shall refer to 
(ini) as "BMn" (the derivative is along the slice normals) 
and to PH)) as "BMt" (the derivative is along the time 
lines) . 

A third slicing condition [l^, 

a^-^lLa{aK-D,p') ee /ii|(ln | det7|)-, (11) 

where Di is the covariant derivative compatible with 
the 3-metric 7^, is also related to BM. We shall call 
it "BMg", as it can be integrated for any hl = 1^l{oi) 
to relate a to the 3-metric determinant. For /i^ — 2 /a, 
BMg integrates to a = f{x) + In | det 7^ |, explaining the 
name "1-1- log slicing". BMn and BMt can be integrated 
only if the shift is zero. 

The geometric specification of BMt and BMg (but not 
BMn) slicing depends on the shift. Here we shall use the 
"fn-driver" 

f3^-(3^(3y = ^iso?{r-m, (12) 
or the "ft-driver" 

/?^=Msa'(r-/^), (13) 
where is the 3-vector defined by 

/, ^ 7^S„- fe - (14) 

in preferred Cartesian coordinates (see Appendix [Cjl . 
With p = 2/3, these are essentially versions of the (im- 
plicit) "Gamma-driver" shift conditions that now domi- 
nate numerical relativity. 

A simple analysis of BMg as a pure gauge system (sim- 
ilar to Appendix El on Minkowski spacetime shows that 
it is well-posed with a fixed shift (see also [l3|), but is 
ill-posed with the fn or ft drivers. We do not consider it 
further. 

In Appendixl^we also show that BMt slicing in combi- 
nation with any shift condition always has both positive 
and negative gauge coordinate speeds. This means that 
on any excision surface of constant radial coordinate r 
there will always be a gauge mode travelling towards in- 
creasing r, and so excision is not possible with this slicing 
condition unless a boundary condition is imposed on the 
gauge at the excision boundary. A similar result holds 
for the ft-driver shift condition. We will mainly use ei- 
ther an algebraic Killing shift (area freezing shift) or the 
fn driver shift. 
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V. COMPATIBILITY OF KILLING 
COORDINATES WITH BMN SLICING 

A. General 

In this section we ask if Killing coordinates exist that 
are compatible with BMn slicing. Although the BMn 
slicing condition is geometrically independent of the shift, 
a{x^,t) only becomes time-independent if the slicing is a 
Killing slicing and the shift is a Killing shift. Substituting 
d = 0, ([!]) and ^ into we find the scalar equation 



(15) 



on the slice. We use the definitions (-L ^)'' =£,"" — (j>n°-, 
V(a^6) = Slid I^ab = ~ -L ^ a^^b to rewrite this equation 
as a partial differential equation for embedding a slice 
with unit normal vector n°: 



(16) 



where -0 = — is related to the gravitational potential 
in a stationary spacetime and 
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ah 



2 I ah 



(17) 



is a symmetric tensor intrinsic to the slice. 

Given that the unit normal vector of a slice t = const, 
is given by 

ua = -aWJ, a = (-VfaiV^t)-i/2, (18) 

the principal part of HI]) is Q°-''^a^bt. As _L°'' is pos- 
itive definite and > 0, two eigenvalues of Q"'' are 
always positive. The third eigenvalue is associated with 
the eigenvector _L^° and is given by 



i? = (ml-i)<^' + V'- 



(19) 



Therefore is elliptic for Z? > and (2+1) hyperbolic 
for D < 0. 

Alcubierre [3l has shown that the BMn slicing condi- 
tion can also be written as the 3-f 1 wave equation 



P'^^VaVfci - 0, 



P'''' = -n'^n" + pLL{a) L''" (20) 



where _L°^, a and n° are as given above. We have 
perturbed this equation around a Killing slicing t, but 
have not been able to identify any lower-order (friction- 
like) terms that would always push 5t locally towards 
^"VaSt ^ or St = 0. We conclude that if BMn shcing 
is really symmetry seeking in some circumstances, as our 
numerical evidence below suggests, this is not because of 
local friction terms, but rather through the mechanism 
by which a solution of the wave equation on a finite do- 
main with a dissipative boundary condition settles to a 
time-independent solution of the Laplace equation. 

The characteristics of the wave equation are null 
surfaces of the "gauge metric" {P~^)ab — —naUij + 
fij^^ J-ab, which is the matrix inverse of P"''. A slice 



evolving under ([20)) can be excised on a boundary ruled 
by trajectories of the Killing vector only if the Killing 
vector is "spacelike" with respect to the gauge metric, 
that is (P~^)afcCC^ > 0- We find that this is once again 
equivalent to < 0. 



B. Schwarzschild spacetime in spherical symmetry 

This subsection reviews and generalises The BMn 
Killing slicing condition in spherical symmetry is 

Pa' = fiLia)a^K (21) 

Using © and ^ to eliminate a and K gives 



R" 13' R' 

' 2— 0, 



R'fiLiR') P R 



(22) 



which has an obvious first integral that can be expressed, 
using © and ([7]), as 



da 



In 



R'- 



c. (23) 



a Ml (a) 

Alternatively, using ([7]) to eliminate (3 from ([^ gives 



R D' 



where 



TV = 2P'2 



3M\ 



R 



D = [^lL{R')~l]R'^ + 1 



2M 



(24) 

(25) 
(26) 



[Here D has the same meaning as in (|T9t .] For given 
^,l{c<) this is a second order ODE for R{r). For the solu- 
tion to be regular for all i? > 0, iV and D have to vanish 
at the same r, which becomes a regular singular point. 
This fixes R and R' at this r, and hence the constant 
c in (P5)) . (|23p can then be solved as a first-order ODE 
for R{r). This means that for any //^(a), there are at 
most a finite number of twice differentiable spherically 
symmetric Killing slicings of Schwarzschild, one for each 
possible regular singular point. The 3-dimensional PDE 
of which (|24p is the reduction to spherical symme- 
try is elliptic for R > Rc and hyperbolic for R < Re- 
in the absence of spherical symmetry, requiring regular- 
ity at the 2-dimensional boundary between elliptic and 
hyperbolic regions would also make the slice more rigid, 
as it does in spherical symmetry. The first integral (j23p . 
however, has no counterpart in the absence of spherical 
symmetry. 

1+log slicing The case of 1-l-log slicing, /i^ = 2/ a 
has been presented in based on earlier work in p^ . 
There are two possibilities for regular singular points. 
One is i? = 2M with R' = 0. This gives a Killing slicing 
where each slice goes through the bifurcation point of 
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the horizon, the lapse is positive in R and negative in L, 
and the shoes never reach P or F. It is not of interest for 
numerical evolutions. 

The other regular singular point is R' — R'^ = ~3 + 
^/lO, R = Rc = M/{AR'^) ~ 1.54057A/. In this solution 
i? — > oo as r — + oo and R Rq from above as r — oo. 
Ro can be found fr'om ([23]) with R' — 0, is given in im- 
plicit form in [1], and is approximately Rq — 1.31241Af. 
Inside the black hole the slices become asymptotically 
tangent to the Killing field and terminate at ij. The in- 
trinsic geometry of each slice becomes a cylinder of radius 
i?o as r -oo (Fig. H]). 

Harmonic slicing Harmonic slicing is the special case 
of BMn slicing with fi^ = I. The regular singular points 
are then R = 2M with either i?' = or it!' = ±1/2. 
The former can be discarded, and the sign in the latter 
is trivial, so that the Killing slices are characterised by 
a = i?' = 1/2 at i? = 2M . These slices stretch from to 
the future singularity R — Q, and so must be used with 
excision. The gauge characteristics are the light cones 
[l7| , so the gauge only requires the excision boundary to 
be spacelike. 

General ijll{oi) Killing slices cannot have an ex- 
tremum of R if they are to be stationary points of some 
slicing condition. From we see that if the Killing 
slices are to approach ij, that is linir^oc i? = > 0, 
the integral 

(27) 

must be finite, for example with /i^ = 2/a. We con- 
jecture that, conversely, if this integral diverges, as with 
fiL = 1, the Killing slices must intersect the future sin- 
gularity. 

Excision and uniqueness One might think that excis- 
ing a BMn Killing slice would make it less rigid, because 
the regular singular point i? = i?c could be excised. This 
is correct if one excises aX R^ < R < 2M and imposes an 
explicit boundary condition on the slicing, for example 
by fixing a at the excision boundary. By function count- 
ing one would expect the value of a at the boundary to 
control the value of the constant c of the slice. However, 
to excise all modes including the lapse gauge modes, the 
excision boundary must be in the region where < 0, 
and so R = Rc must be on the slice. The only possi- 
ble Killing endstate of the slicing is then the unique one 
derived above. 

VI. VACUUM BLACK HOLE EVOLUTIONS 
A. Method 

To see empirically if generic black hole evolutions are 
attracted to the Killing states we have characterised 
above, we have carried out numerical evolutions of the 
Schwarzschild spacetime in spherical symmetry, using 
BMn 1+log slicing. 



We can take advantage of the fact that this metric is 
known in closed form to evolve only the coordinates on 
the known spacetime, see Appendix|Xl There is no global 
coordinate system that covers wormhole slices and which 
is also Killing. Therefore, in pure gauge evolutions of 
wormhole slices stretching from to i^, we restrict to 
slices with a discrete "left-right" isometry through the 
coordinate sphere r = 0, so that we only evolve explicitly 
on F and R, where KS coordinates can be used, with a 
boundary condition at r = representing the isometry. 

Even this does not work for slices which go through the 
horizon bifurcation 2-sphere (where KS time and similar 
Killing time coordinates are — oo), and so for such slices 
we need to evolve the Einstein equations in the NOR for- 
mulation, see Appendix |B1 In all other cases, plots are 
from pure gauge evolutions, but we have verified that 
our results are replicated in evolutions of the full Ein- 
stein equations in the NOR formulation. The evolutions 
described here all use the fn shift condition (fT2)) except 
otherwise stated. 



B. With excision boundary 

As initial data for the geometry and the coordinates 
we have considered: 

la) KS sfice, KS lapse, KS shift, area radius; 
lb) KS slice, KS lapse, zero shift, area radius; 

2) A closed form asymptotically cylindrical slice, unit 
lapse, zero shift, area radius, see Appendix lA 21 

3) The Hannam slice, lapse, shift, all in area radius, 
see Appendix I A 21 

We first evolved with area locking (that is. Killing) 
shift (which is determined algebraically so that the ini- 
tial value of the shift listed above is irrelevant) and exci- 
sion. We excised at i? = 1.54M, which is just inside the 
maximal excision radius R = Rc — 1.54057M for which 
all modes are outgoing. We find that 1) and 2) approach 
the Killing state, and 3) remains there. This is demon- 
strated in Fig. [51 and indicates that the Killing state has 
a significant basin of attraction. 

When combined with the fn shift driver, in la) the co- 
ordinates r are pushed out of the black hole and further, 
lb), and 2) again settle down to the Hannam endstate, 
and 3) remains there. With lb), the excision radius ini- 
tially has to be i? ~ 1.3Af or the excision surface at con- 
stant r will be pushed out so far before it reaches steady 
state that there a gauge mode is ingoing. 

C. With isometry boundary condition 

We begin with NOR evolutions starting from the time- 
symmetric wormhole slice that goes through the bifurca- 
tion 2-sphere R = 2M. 

We first use spatial coordinates in which i'^ is repre- 
sented by the point r = (the "puncture"), see Ap- 
pendix [B| Evolutions with this method reproduce the 
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FIG. 5: The distance of the lapse from the Kilhng endstate 
over the range from the excision boundary R = 1.54 (just 
inside the regular singular point) out to _R = 21.54, with area 
locking shift. The power law decay indicates ||a — ckKiiiing|| ~ 
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FIG. 6: The distance of the lapse from the Kilhng endstate 
from the excision boundary with approximately the same lim- 
its as in Fig. [S] using the fn shift driver (|12p . The power law 
decay again indicates ||q — ciKiiiing|| t^^ ■ 



behaviour described in numerical error changes the 
topology and the evolution settles down to the asymp- 
totically cylindrical Killing state. 

If we evolve the same initial data in spatial coordinates 
that resolve the wormhole (see Appendix [B|, we see the 
slices begin to form a cylinder at radius i?o , but at reason- 
able resolution constraint violation in the Einstein code 
makes the result unreliable soon after. 

Pure gauge evolutions with wormhole initial data that 
lie to the future of the bifurcation 2-sphere (so that 
R < 2M at the throat) and a discrete isometry boundary 
as described in Appendix [Xl are more stable. At the isom- 
etry boundary (where R is minimal) the lapse quickly col- 
lapses. In low resolution evolutions, the lapse collapses 
starting at the minimal R (at the isometry boundary), 



and a cylinder of radius i?o forms with proper length 
increasing linearly in time (Fig. [7]). 

D. Gauge shocks 

However, higher resolution (for example Ar = M/50) 
evolutions show that low resolution only hides the for- 
mation of a gauge shock where K forms a large nega- 
tive peak and a" forms a positive peak, at i? ~ 1.5M. 
This does not seem to happen exactly at Rc (we varied 
fiL{(^) to check this), and so we do not think that it is 
a kink instability related to the regular singular point of 
the Hannam slice. Neither is there any indication that 
the slice has become null. An ODE mechanism by which 
K < makes a grow is also ruled out as not all initial 
data where K < shock. 

Rather, we think we see a gauge shock of the type 
described by Alcubierre [l^, [l^l- Note that the lapse 
speeds expressed in terms of proper distance / per coor- 
dinate time t, relative to the time lines, are —f3±a^^L = 
— /3±\/2a, so that a gauge wave propagating "left", from 
high to low a is expected to steepen. By contrast, the 
wave propagating "right" and forming the cylinder ap- 
pears to be stable and translating with constant speed 
dl/dt without changing its shape much. Alcubierre notes 
that for the particular choice /iL = 1 + k/a^ with k > 
the pure gauge system is linearly degenerate, and we have 
tried this fiL, but shocks still form, also in agreement with 
Alcubierre's numerical observations. Alcubierre argues 
that gauge shocks are generic for evolved gauge condi- 
tions. 

Although the NOR evolutions of the time-symmetric 
slice are less reliable, they suggest that evolutions shock 
when a has a local minimum not at the isometry bound- 
ary (Fig. [9]). They also suggest that with a — 1 ini- 
tially the slicing never shocks (Fig. [5]) . This agrees with 
the standard numerical literature where the puncture 
data are approximately the time-symmetric slice through 
Schwarzschild and the initial lapse is one. It also agrees 
with the evolution by Brown [l9| of these particular ini- 
tial data. It seem plausible that initial data in a neigh- 
bourhood also do not develop shocks, but we have not 
investigated this. 

We note that the shift remains regular during the blow- 
up, and the same qualitative picture occurs with proper 
distance radius, zero shift, or fn driver shift. 

With BMt slicing and the ft (not fn) shift driver, we 
see the same gauge shock in both NOR and pure gauge 
evolutions, but it seems to form earlier and even at low 
resolution, so that we never see formation of a cylinder 
before the code crashes. 



VII. SCALAR FIELD COLLAPSE EVOLUTIONS 

We now consider the behaviour of collapse simulations 
with BMn 1-flog slicing. As a toy model we consider 
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time = 1 .20000 




FIG. 7; Snapshots of R and a against proper distance r from 
an evolution of an isometric slice. The throat of the slice (ini- 
tially at i? = 1.5M) is gradually stretched so that it becomes 
an infinitely long cylinder. The radius of the cylinder agrees 
with that computed in Q. Note that low numerical resolu- 
tion effectively smears out a gauge shock travelling left, so 
that this is not a correct continuum solution. 



time = 83.75000 




8 10 



FIG. 8: The K = time symmetric slice through the bi- 
furcation surface of Schwarzschild, evolved with BMn l-|-log 
slicing, with a — 1 initially. We show a snapshot of K and 
a against proper distance radius. The edge at r ~ 7 in this 
snapshot moves to the right, and leaves behind a cylinder of 
constant R and K with a ~ 0. 



a spherical scalar field. We impose spherical symmetry 
and use proper distance as the radial coordinate. The 
metric thus takes the form (j4|) with 7 = 1. Details of the 
numerical implementation and the initial data are given 
in Appendix [P] 

The initial data are chosen as a moment of time sym- 
metry. The scalar field separates into an ingoing pulse 
and an outgoing pulse. With the chosen parameters, the 
ingoing pulse collapses to form a black hole, with an ap- 
parent horizon first forming at t — 6.9. The final mass of 
the black hole is 1.0. Figs. [TUl and [TT] show respectively 
a and K at t = 14 and in the range < r < 3. Note 



FIG. 9: As in Fig. |8l but with a not constant on the initial 
slice. The wave on the left travels left and is steepening, about 
to form a gauge shock, with large negative K. The wave on 
the right travels right and is also steepening: note a" > 
there. 




FIG. 10: Plot of the lapse a against proper distance r at 
t = 14 in scalar field collapse without excision. 



the sharp features in both these quantities near r — 1.5. 
These features become ever sharper and cause the code 
to crash not long after the time of these graphs. This 
pathology is again a gauge shock. Neglecting the shift 
the principal part of the the evolution equation for K is 
K = —a". Combining this with the BMn equation yields 
a nonlinear wave equation for the lapse whose principal 
part is a = 2aa" . The modes of this equation travel with 
speeds i^/2a. Thus if one has an inner region where the 
lapse has collapsed, then left-moving gauge waves pile up 
on the boundary of this region, giving rise to a shock wave 
in a which in turn (through the equation K = —a") will 
induce a shock wave in K. Fig. [T^] shows a" a,t t — 14 
and in the range < r < 3. Note that this quantity also 
has a sharp feature near r = 1.5. 

In this simulation, the pathological behaviour is inside 
the horizon. This suggests that we might be able to avoid 
the pathology by using excision. Figs. [HI and [T51 respec- 
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FIG. 11; Plot of K against proper distance r at t = 14 without 
excision. 




0.5 1 1.5 2 2.5 3 



FIG. 12: Plot of a" against distance r at t = 14 without 
excision. 

lively present the values of a and K for a simulation done 
using excision. Here the simulation is run until t — 60.5 
and all quantities are plotted as functions of the area ra- 
dius R rather than proper distance r. It is excision that 
allows the simulation to be run this long because the ex- 
cised grid contains no regions of negative K which caused 
the non-excision simulation to crash. Furthermore, this 
late in the simulation, these quantities have asymptoted 
to the static values described in Q . This is illustrated in 
Fig. [Til which contains two plots of a as a function of R. 
The solid line is the a given by the endstate of the exci- 
sion collapse simulation, while the dotted line represents 
the Killing lapse given by an integration of the ODEs of 

i- 

VIII. CONCLUSIONS 

Killing endstates We have explained why it is possi- 
ble in evolutions of black holes that all metric coefficients 
become time-independent without either slice stretching 
or collapse of the lapse. We have reviewed the Bona- 
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FIG. 13: Plot of K against area radius R ai t = 60.5 with 
excision. 




5 10 15 ^ 25 30 35 40 

FIG. 14: Plot of Q against area radius R a.t t = 60.5 with 
excision (solid line) and the exact Killing lapse (dashed line). 



Masso slicing conditions, and have derived a mixed ellip- 
tic/hyperbolic PDE on the slice that characterises Killing 
endstates of the BMn family of slicing conditions. Nu- 
merically, we have shown that spherical BMn slicings of 
the Schwarzschild spacetime are attracted to the Killing 
endstate from nearby initial data. We do not fully under- 
stand the mechanism for this. Initial data further away 
also appear to be attracted to the Killing state, but on 
closer inspection this is true only at low numerical reso- 
lution. 

Gauge shocks Increasing the resolution reveals that in 
the continuum the 1+log BMn slicing generically devel- 
ops gauge shocks of the type described by Alcubierre [l^ , 
where the speed of gauge waves associated with the slic- 
ing increases with the lapse, so that gauge waves moving 
from large to small lapse steepen. The only initial data 
set we have examined that does not form a gauge shock 
with BMn 1-l-log slicing is the time-symmetric wormhole 
slice through Schwarzschild spacetime with unit initial 
lapse, although we suspect that there is at least a neigh- 
bourhood of such data. More numerical work is required 
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to explore this. 

It may be that gauge shocks would also occur in bi- 
nary black hole simulations with l-|-log BMn slicing in 
the continuum, but that they are suppressed by low res- 
olution inside the black holes. By contrast, in collapse 
simulations the central region is typically adequately re- 
solved, and in fact recent work where the collapsing re- 
gion is never excised seems to require large dissipation 
for stability 

Excision We find that in both collapse and vacuum 
simulations the gauge shock can typically be avoided by 
excising just inside the apparent horizon. There seems 
to be no clear awareness in the literature that such a 
boundary still has an incoming gauge mode, and that 
the resulting continuum problem is ill-posed. 

Confirming this, our vacuum (both Einstein and pure 
gauge) evolutions in spherical symmetry do not converge 
and often blow up when an incoming mode at the excision 
boundary is neglected. By contrast, our collapse code, 
which uses different numerical methods, does not seem to 
mind. There is also an explicit claim that 3D binary black 
hole evolutions converge equally well with and without 
excision [lollTl|. 

We have also re-derived the previously known [l3| fact 
that if gauge drivers are not of the form a + (3^a^i = . . . 
and /?* + (3^ P'^j — ■ ■ ■, full excision is not possible at any 
radius. 

Nature of the Killing endstate of BMn slicing If there 
is to be no incoming gauge mode at an excision boundary, 
the equation obeyed by the Killing slice has a transition 
from elliptic to hyperbolic. Requiring regularity there 
makes the slice more rigid, and in spherical symmetry 
makes it unique. The same is true if the slice has no 
excision boundary. 

We have clarified that this unique Killing BMn slicing 
of a Schwarzschild black connects spacelike infinity out- 
side the black hole to future timelike infinity inside the 
black hole, where it becomes asymptotically cylindrical. 
As pointed out independently by Brown ^IS] , initial data 
which connect two asymptotically flat regions through a 
wormhole cannot evolve to this endstate in the contin- 
uum, although numerical under-rcsolution gives the false 
impression that the topology jumps. In the continuum 
evolution, the wormhole stretches into a cylinder whose 
length grows without bound. 

Comments on 3D evolutions In our investigation, we 
have identified three problems with current gauge choices 
in 3D numerical evolutions of collapse and black holes 
with a currently favoured slicing condition, BMn l-|-log 
slicing: 1) wormhole data do not admit a BMn Killing 
endstate; 2) excision close inside the apparent horizon re- 
quires explicit boundary conditions for the gauge; and 3) 
coordinate shocks form generically. None of these prob- 
lems have been noted in the binary black hole literature, 
but we believe that this is only because of limited resolu- 
tion, and that they will become apparent as a failure of 
convergence or instabilities at sufficiently high resolution. 
There are, however, simple ways around these problems: 



• Wormhole initial data for eternal black holes ending 
at «° should be replaced by initial data that are 
asymptotically cylindrical and end at ij. 

• Continuum boundary conditions should be imposed 
explicitly at excision boundaries for any incoming 
gauge modes. 

• The initial lapse should be chosen such that gauge 
shocks do not form. This will require more em- 
pirical studies in 3D. In collapse without excision, 
changing to a smoothly collapsed lapse profile once 
an apparent horizon has formed may be helpful. 

Final remarks Finally, two technical developments 
given in the appendix may also be of interest to the 3D 
community. By characterising pure gauge as the evolu- 
tion of a coordinate system on a background spacetime 
given as = X^{x^,t)^ we have been able to check 
strong hyperbolicity of the gauge and calculate the gauge 
speeds without reference to a formulation of the Einstein 
equations. By re-defining the vector auxiliary variable 
of the NOR and BSSN formulations (following HJ), we 
have made them easier to use with non- Cartesian coor- 
dinates or multiple coordinate patches. 
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APPENDIX A: PURE GAUGE EVOLUTIONS 

1. General equations 

Prescriptions for the lapse and shift can be tested with- 
out the specific stability problems associated with any 
formulation of the Einstein equations by evolving coordi- 
nates (ijX') on a spacetime given a priori in coordinates 
with metric g^,j. The ADM evolution equations are 
replaced by the definition ([T]) of d/dt in terms of the lapse 
and the shift, 

= an^' + P'X^^, (Al) 

where are the X^ components of the unit normal to 
the t-slices. The 3-metric is given by 

(A2) 

and = Cp^ij - 2aKtj gives 

K,, = X''^,X9>.. + l{x''.rX\X 

+X\,X\jn^ - X^^,X\y^g^,^^, (A3) 



11 



where the unit normal is given, up to normaHsation, 
by 



(A4) 



In spherical symmetry is defined by (d/dr)'^na ~ 0, 
and in preferred coordinates (T, K) 



R' T' 
nr = riR^—. 



Initial data 



(A5) 



We use KS coordinates for the background 
Schwarzschild spacetime. Simple closed form slices 
which become asymptotically cylindrical at R = Rq may 
be constructed by making the ansatz T = t + F{R) with 



F\R] 



-LiR-Ro) 



(A6) 



Constructing 1+log Killing data is straightforward in 
area gauge where r — R. We note that proper distance 
I is dl — ^/JdR in area gauge r — R, while R' = dR/dl 
in proper distance gauge r — I. Therefore, we obtain 
an algebraic expression for 7 in area gauge by replacing 
R' by \/7 in (|23p . A single numerical integral must be 
performed to compute F(r). 

A surface with the required isometry in the F re- 
gion, with a radial coordinate r in which the isometry 
is r — !■ — r, can be given in Schwarzschild coordinates 
(77, i?) as R{—r) — R{r) and r]{—r) - rj^ = —{ri{r) — rj^). 
Here 77 = 77* is the reflection surface. Note it is time- 
like in F. The last condition must be translated into a 
relation between T(— r) and T(r) using rj — T — 4>{R) 
where (f>{R) = 2Af ln(l — 2M/R). Furthermore, we want 
the Kerr-Schild time T to be smooth through R — 2M, 
while 4'{R) is not. An ansatz with these properties is 



T(r) = C{r)r + [1 + D{r)](l){R{r)), 



(A7) 



where C(r) and D{r) are smooth odd functions and D{r) 
is at r = and —1 on and outside the horizon. To 
preserve the isometry and the adapted radial coordinate 
in the evolution the lapse and shift must obey a(— r) — 

air), Pi-r) = -/3(r). 



3. Characteristic analysis 

The spherical gauge evolution system is not quasilin- 
ear, because the equations for R and T are nonlinear in 
R' and T'. Hence in order to analyse the hyperbolic- 
ity of the system one must explicitly linearise, and apply 
weights to equations and variables [2^ . The principal 
part of the evolution of the coordinates is then 



ST 



~ i 9RR - [n^f \ 5R' + T75a + T'5f3{M) 



SR ~ lg^^ + ^n''f \6T' 



V7 



\SR' 



by 



+n"Sa + T'S(3. (A9) 
For BMn slicing with ft shift driver, this is completed 



Sa ■ 
S(3- 



7 



7 



[urST" + utSR") + 136a', (AlO) 
[n'^dT" + n^6R"). 



3/2 



(All) 



The system is diagonalisable with real characteristic 
speeds dr/dt 



^/7 



-/3± 



1 (^^±^^2 + wsrj, 



(A12) 
(A13) 



and so is strongly hyperbolic. 

The principal part of BMn with the fn shift driver 
is the same except that 5/3 = ... -f pSfi' . This system is 
also strongly hyperbolic with characteristic speeds 



(AM) 
(A15) 



For BMn with area locking shift the principal part be- 
comes 



5T 



f35T' - —Sa, 

Tlx 



1 



ST" + (3Sa'. 



(A16) 
(A17) 



In this case the system is also strongly hyperbolic and 
has characteristic speeds (|A12p . 

The BMt slicing condition (fTU)) gives rise to strongly 
hyperbolic pure gauge systems with speeds 



(AI8) 



with either area locking shift, ft shift driver [adds speeds 
(|A13|) ] or fn shift driver [adds speeds (|A15|) ] . 

Constructing boundary conditions for a system which 
is not quasilinear is a difficult task which will not be 
discussed here. Numerically, we have simply frozen all 
fields at the outer boundary, but we have moved the outer 
boundary so far out that in the continuum it does not 
affect the results shown here. 
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APPENDIX B: WORMHOLE INITIAL DATA FOR 
THE EINSTEIN EQUATIONS 

To create a numerical evolution in spherical symmetry 
that is similar to the BSSN moving puncture evolutions 
in 3D, we use the well-known isotropic radial coordinate 



(Bl) 



with range Q < r < oo. We use the correct reduction 
to spherical symmetry of the NOR system described in 
Appendix [Cl stagger the grid around r = 0, and im- 
pose as boundary conditions at r = that 7rr, R and 
a are even in r, and and /?'' are odd. These condi- 
tions hold for the puncture initial data, and it is easy 
to see that they are compatible with the time evolution. 
If r = is a regular centre of spherical symmetry, then 
(— r, 9, ip) represents the same point as (r, tt — 0, + tt), 
and these conditions follow from spherical symmetry and 
regularity. If r = represents then r < is simply 
not part of the spacetime, the even/odd conditions do 
not follow from spherical symmetry, and their meaning 
is unclear. However, for finite differencing purposes they 
are equivalent to finite differencing across the puncture 
X — y = z — Q in 3D as if ii was a regular point, which 
is what is done in 3D puncture evolutions. 

In order to resolve both sides of the wormhole, we use 
the symmetric radial coordinate 



R(r) = v/r2 +4M2, 



(B2) 



with range — cx3 < r < oo, and impose isomctry boundary 
conditions at r = 0. 



APPENDIX C: REDUCTION OF THE NOR 
FORMULATIONS TO SPHERICAL SYMMETRY 

To our knowledge, little numerical work using the 
BSSN or NOR formulations in spherical symmetry has 
been published, and therefore we would like to point out 
a technical detail in the reduction to spherical symmetry. 
fi can be defined as a true 3-vector by introducing a flat 
auxiliary connection [2l| , so that 



fa = l^^'^dab 



(Cl) 



Here V is defined so that its connection coefficients van- 
ish in preferred Cartesian coordinates x*. In an evolution 
with a single Cartesian coordinate patch these are the 
only coordinates, and Vq reduces to the partial deriva- 
tive, but if multiple Cartesian patches are used, one of 
them is preferred. 

In spherical symmetry, Va is defined as the covariant 
derivative compatible with the fiat metric ds^ — dr^ + 
r^dil'^. The result is 

f^ = 7k + l(7lL^l]^P(2k + 22^] (C2) 

Jrr r ) 2 y^rr IT J 



where 7t = /r"^ ^ we have written explicitly instead 
of just 7 as elsewhere in this paper, and /e = = 
because of spherical symmetry. Local regularity at the 
origin r — Q requires ^rr — It = 0{r^) and fr — 0{r), 
and so fa is a regular vector field. [Naively calculating 
/,■ using partial derivatives in (r, 6, ip) gives a different, 
singular result.] 



APPENDIX D: SPHERICAL EINSTEIN-SCALAR 
CODE 

We choose the collapsing matter to be a massless, min- 
imally coupled scalar field ip so the equations of motion 
become 



Rab = Vai'^bi^, 



(Dl) 
(D2) 



We will find it helpful to introduce the quantities A, P, 
H, s and w by 

A = Krr- \K, (D3) 

P = i,\ (D4) 

H - n"Va^, (D5) 

s = a', (D6) 

w = R', (D7) 

and also to define the quantities a and s by 

a = lna, (D8) 

s = s/a. (D9) 

The spatial metric evolves by 

dtlij = -~2aKij + Cpjij. (DIO) 

Then noting that jrr = 1 we find that the two compo- 
nents of (jPToli become 



= -a{A+ ^K) + (3', 
R = l3w + aR{\A-\K). 



(Dll) 
(D12) 



The first of these equations can immediately be inte- 
grated to yield 

/3 = J a{A+ ^K) dr. (D13) 

Differentiating the definition of P with respect to time 
we obtain 

P ^ (3P' + sU + a[li' + P{lK + A)] , (D14) 
while the wave equation for yields 

U^pn' + sP + aiKU + P' + . (D15) 
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The momentum and Hamiltonian constraints yield re- 
spectively 

A' - - nP + IK', (D16) 

^' = ^ - I (tf + + lA- - IK') .(D17) 

The momentum constraint enables us to differentiate 
()D12p with respect to r yielding 

Pw' + sR{lA-\K)^laRPIV. (D18) 

The BMn slicing condition is 

a = fi~s- 2K, (D19) 

which when differentiated with respect to r yields 

S = + aS {A + ^K) - 2K'. (D20) 

The extrinsic curvature evolves by 

£tK\ = CpK^b - D^Dba 
+ai^^^R''b + KK'^b- D^'ijDbiP), (D21) 

where Da and *^'^^P''h are respectively the derivative op- 
erator and Ricci tensor associated with the spatial metric 
Jab- Eq- (|D2ip yields the evolution equations 

k^f3K'-s'-^ + a {\K^ + \A^ + n2|D22) 

+a (^itf - ip2 + \A^ - + KA+ !^^yB2S) 

We evolve this system in one of two different ways de- 
pending on whether the numerical inner boundary is a 
regular centre or an excision boundary. 

If the inner boundary is a regular centre then the vari- 
ables to be evolved are {R, P,Il,w,a,s, K) which are 
evolved using Eqs. (|D12ID14ID15ID18ID19ID20llD22)) re- 
spectively. In the evolution equations, the quantities a, s 
and s' are regarded as derived from the evolved variables 



a and s. The shift /3 is determined from (|D13p and the 
condition that /3 vanish at the centre. Similarly A is com- 
puted by integrating (|D16p outward from the centre and 
using the fact that A vanishes at the centre. Boundary 
conditions at the centre simply follow from smoothness 
of the metric which requires that _R, P and s vanish there, 
that the spatial derivatives of 11, a and K vanish, and 
that w = 1. 

If the numerical inner boundary is an excision bound- 
ary, then the variable A is not computed by the momen- 
tum constraint but is instead evolved using (|D23[) . The 
integration constant for ID13|) is determined by having 
the time derivative of the area radius R vanish at the 
apparent horizon. At the excision boundary, all modes 
are outgoing, except for a single incoming gauge mode 
s — ^/2aK. As a boundary condition, its value is pre- 
scribed as the value it had at the previous time step. All 
other modes are evolved at the excision boundary point; 
however spatial derivatives at that point are computed 
using one sided differences rather than the centred dif- 
ferences used at the other points [1^. 

For initial data we choose a moment of time symmetry 
so that K and 11 vanish. The initial value of P is given 
by 

P = ar^ exp(-rVcr^) (D24) 

where a and a are constants, and the initial value of R is 
solved for by using the Hamiltonian constraint. The lapse 
a is given an initial value of 1 . The grid points are equally 
spaced in r. Spatial derivatives are evaluated using stan- 
dard second order centred differences and time evolu- 
tion is done by the iterative Crank-Nicholson method. 
A marginally outer trapped surface occurs where the 
derivative of R along the outgoing null direction vanishes. 
In terms of the variables used here, that corresponds to 
the vanishing of w + R{^A ~ ^K). If the apparent hori- 
zon (the outermost marginally outer trapped surface) oc- 
curs at gridpoint i (and if we choose to excise), then we 
place the excision boundary at gridpoint 9i/10. The con- 
stants determining the initial value of P were chosen to 
be a = 0.05 and a' = 8. The number of spatial grid- 
points was chosen to be 6401 and the initial range of r 
was chosen to be < r < 200. 
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